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THE  APEX  METHOD  IN  IMAGE  SHARPENING  AND  THE  USE  OF 
LOW  EXPONENT  LEVY  STABLE  LAWS 


ALFRED  S.  CARASSO* 

Abstract.  The  APEX  method  is  an  F FT- based  direct  blind  deconvolution  technique  that  can 
process  complex  high  resolution  imagery  in  a few  minutes  of  cpu  time  on  current  desktop  platforms. 
The  method  is  predicated  on  a restricted  class  of  shift-invariant  blurs  that  can  be  expressed  as 
finite  convolution  products  of  two-dimensional  radially  symmetric  Levy  stable  probability  density 
functions.  This  class  generalizes  Gaussian  and  Lorentzian  densities  but  excludes  defocus  and  motion 
blurs.  Not  all  images  can  be  enhanced  with  the  APEX  method.  However,  it  is  shown  that  the  method 
can  be  usefully  applied  to  a wide  variety  of  real  blurred  linages,  including  astronomical,  Landsat 
and  aerial  images,  MRI  and  PET  brain  scans,  and  scanning  electron  microscope  images.  APEX 
processing  of  these  images  enhances  contrast  and  sharpens  structural  detail,  leading  to  very  noticeable 
improvements  in  visual  quality.  The  discussion  includes  a documented  example  of  non  uniqueness 
where  distinct  point  spread  functions  produce  high-quality  restorations  of  the  same  blurred  image. 
Significantly,  low  exponent  Levy  point  spread  functions  were  detected  and  used  in  all  the  above 
examples.  Such  low  exponents  are  exceptional  in  physical  applications  where  symmetric  stable  laws 
appear.  In  the  present  case,  the  physical  meaning  of  t hese  Levy  exponents  is  uncertain. 

Key  words,  image  deblurring;  blind  deconvolution;  direct  methods;  electronic  imaging  systems; 
low  exponent  stable  laws:  APEX  method;  SECB  method;  non  uniqueness;  astronomical,  Landsat, 
and  SEM  images;  MRI  and  PET  brain  scans. 

AMS  subject  classifications.  35R25,  35B60,  60E07,  68U10. 


1.  Introduction.  The  APEX  method  is  an  FFT-based  direct  blind  deconvolu- 
tion technique  introduced  by  tire  author  in  [9].  The  significance  of  the  present  paper 
lies  in  the  successful  use  of  that  method  in  sharpening  a wide  variety  of  real  blurred 
images , as  opposed  to  the  synthetically  blurred  images  discussed  in  [9].  The  rea- 
sons behind  these  successful  applications  are  not  fully  understood.  Not  all  images 
can  be  usefully  enhanced  with  the  APEX  method.  The  present  paper  is  essentially 
self-contained  and  may  be  read  independently  of  [9]. 

Blind  deconvolution  seeks  to  deblur  an  image  without  knowing  t he  point  spread 
function  (psf)  describing  the  blur.  Most  approaches  to  that  problem  are  iterative 
in  nature.  Because  non  uniqueness  is  compounded  with  discontinuous  dependence 
on  data,  such  iterative  procedures  are  not  always  well-behaved.  When  the  iterative 
process  is  stable,  several  thousand  iterations  may  be  necessary  to  achieve  useful  re- 
constructions. However,  as  shown  in  [9],  by  limiting  the  class  of  blurs,  non  iterative 
direct  procedures  can  be  devised  that  accomplish  blind  deconvolution  of  512  x 512 
images  in  a few  minutes  on  current  desktop  platforms. 

The  APEX  method  assumes  the  image  g(x,y)  to  have  been  blurred  by  a restricted 
type  of  shift-invariant  psf  h(x,  t /),  one  that  can  be  expressed  as  a finite  convolution 
product  of  2-D  radially  symmetric  Levy  stable  probability  density  functions.  Such  so- 
called  class  G psfs  include  Gaussians,  Lorentzians,  and  their  convolutions.  However, 
the  class  G also  excludes  defocus  and  motion  blurs,  and  convolutions  of  such  blurs 
with  Gaussians  and  Lorentzians. 

The  synthetically  blurred  images  g{x,,y)  used  in  [9]  were  created  by  numerical 
convolution  of  sharp  images  f(x,y)  with  class  G psfs  h.(x,  y).  Such  blurred  images 
necessarily  obey  the  convolutional  model  g(x,  y)  = h(x,  y)<  > f(x,y)  + noi.se,  on  which 
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the  APEX  method  is  predicated.  In  a real  image,  the  blur  need  not  be  radially  sym- 
metric nor  shift-invariant,  and  may  otherwise  be  poorly  approximated  by  an  element 
of  G.  More  fundamentally,  the  blurring  operator  may  not  even  be  linear.  Applicabil- 
ity of  the  APEX  method  to  a given  real  image  is  far  from  obvious.  Therefore,  useful 
sharpening  of  any  such  image  with  an  APEX-detected  psf  is  always  instructive. 

Stable  distributions  are  the  natural  generalization  of  the  Gaussian  distribution. 
Their  theory  was  developed  by  Paul  Levy  in  the  1930’s  in  connection  with  his  work  on 
the  Central  Limit  Theorem,  [13].  In  the  simplest  radially  symmetric  case,  these  dis- 
tributions are  characterized  by  an  exponent  /3,  0 < (3  < 1.  with  (3  = 1 corresponding 
to  the  Gaussian  distribution,  and  [3  = 1/2  corresponding  to  the  Cauchy  or  Lorentzian 
distribution.  Because  stable  distributions  have  infinite  variance  when  (3  < 1,  their  ap- 
pearance in  physical  contexts  sometimes  poses  philosophical  difficulties.  In  the  present 
case,  use  of  such  psfs  as  the  framework  for  the  APEX  method  is  motivated  by  the 
important  role  Levy  densities  appear  to  play  in  numerous  imaging  systems.  This  is 
documented  in  section  2.  When  the  APEX  method  is  applied  to  a given  image  in  the 
manner  described  below,  a Levy  psf  with  a specific  value  of  (3  is  necessarily  detected. 
That  value  of  ,3  may  not  be  indicative  of  the  actual  physical  process  that  created  the 
image.  This  is  true  even  if  deblurring  with  the  detected  psf  significantly  improves  the 
image.  As  shown  in  section  4,  there  are  in  general  infinitely  many  distinct  values  of  (3 
that  can  produce  useful  reconstructions  from  the  same  blurred  image.  In  some  cases, 
the  usefully  enhanced  image  may  not  have  been  blurred  by  a class  G psf  to  begin 
with.  In  other  cases,  APEX  processing  does  not  significantly  improve  the  image. 

Below,  we  exhibit  ten  images  where  APEX  processing  provided  noticeable  im- 
provement. These  examples  encompass  such  diverse  imaging  applications  as  astro- 
nomical, Landsat,  and  aerial  images,  MRI  and  PET  brain  scans,  a scanning  electron 
microscope  image,  a face  image,  and  other  types  of  interesting  images.  In  some  cases, 
the  improvement  is  due  primarily  to  an  increase  in  contrast.  In  other  cases,  there  is 
demonstrable  sharpening  of  structural  detail  in  addition  to  increased  contrast.  In  all 
cases,  the  change  in  image  quality  is  more  than  cosmetic,  as  APEX  processing  signif- 
icantly alters  the  image  Fourier  transform.  It  is  noteworthy  that  low  exponent  stable 
laws,  with  / 3 1/2,  were  detected  and  used  to  deblur  all  of  the  images  shown  below. 
Such  /3-values  are  exceptional  in  physical  contexts  where  radially  symmetric  Levy 
densities  appear.  Whether  or  not  these  values  have  a physical  origin  cannot  be  as- 
certained in  the  present  case.  Moreover,  the  APEX  detection  procedure  may  not  be 
well-founded.  Nevertheless,  the  fact  remains  that  the  use  of  such  psfs  produced  valu- 
able restoration  of  real  imagery  from  important  fields  of  science  and  technology.  To 
tin1  author’s  knowledge,  this  application  of  sub-Cauchy  stable  laws  in  image  processing 
is  new  and  unanticipated. 

In  recent  years,  there  has  been  considerable  interest  in  image  processing  tech- 
niques that  can  be  formulated  as  initial  value  problems  in  nonlinear  partial  differ- 
ential equations.  An  instructive  survey  of  these  developments  may  be  found  in  [10]. 
In  particular,  several  novel  approaches  to  image  enhancement  have  been  devised, 
based  on  integrating  well-posed  anisotropic  nonlinear  diffusion  equations.  In  contrast, 
the  APEX  method  centers  around  ill-posed  continuation  in  linear  fractional  diffusion 
equations.  The  results  obtained  here  appear  to  compare  favorably  with  what  is  fea- 
sible with  nonlinear  methods,  and  they  indicate  the  APEX  method  to  be  a useful 
addition  to  this  developing  methodology. 

2.  Imaging  systems,  Levy  processes,  and  the  class  G.  The  occurrence  and 
analysis  of  Levy  processes  in  the  physical  sciences  are  subjects  of  significant  current 
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interest.  See  [1],  [2],  [4],  [27],  [29],  [30],  [33],  and  references  therein.  An  important 
special  case  involves  2D  radially  symmetric  Levy  stable  densities  h(x,y),  implicitly 
defined  in  terms  of  their  Fourier  transforms  by 

(1)  //(£,//)  = / h(x,y)e~2iri{ix+riy)dxdy  = e-Q(C2+^)d,  a > Q.  o <f3<\. 

The  cases  (3  — 1 and  (3  = 1/2  correspond  to  Gaussian  and  Lorentzian  (or  Cauchy) 
densities  respectively.  For  other  values  of  /?,  h(x,y)  in  (1)  is  not  known  in  closed  form. 
When  f3  — 1,  h,(x,y)  has  slim  tails  and  finite  variance.  For  0 < (3  < 1,  h.(x,y)  has  fat 
tails  and  infinite  variance.  As  noted  in  [33],  there  are  examples  in  science  where  the 
occurrence  of  a stable  law  can  be  deduced  from  ‘first  principles’  in  terms  of  physical 
mechanisms  that  do  not  involve  the  parameter  / 3 explicitly.  One  such  instance  is  the 
Holtsmark  distribution  describing  the  gravitational  field  of  stars  [13].  There,  mathe- 
matical analysis  reveals  the  value  (3  = 3/4.  Such  cases  must  be  distinguished  from  the 
many  other  cases  where  empirically  obtained  data  with  fat  tails  are  fitted  to  a Levy 
law,  and  the  exponent  f3  is  inferred  from  these  data.  Given  the  limitations  of  physical 
measurements,  such  empirically  established  Levy  processes  do  not  have  the  degree  of 
scientific  legitimacy  that  attends  the  Holtsmark  distribution.  The  considerations  of 
the  present  paper  generally  lie  in  this  weaker  scientific  realm.  Nevertheless,  as  will  be 
seen  below,  techniques  derived  from  such  considerations  turn  out  to  be  effective. 

Image  intensifiers,  CCDs,  and  numerous  other  electronic  devices  are  used  in  a 
wide  variety  of  astronomical,  industrial,  biomedical,  military,  and  surveillance  imaging 
systems.  See  [3],  [11],  [12],  [14],  [20].  Each  such  device  has  a psf  h(x,  y),  characterizing 
that  device’s  imaging  properties.  The  psf  is  a probability  density  function  since  it  is 
non-negative  and  integrates  to  unity.  Use  of  such  a device  to  image  an  object  f(x,  y) 
produces  a blurred  image  g(x,y)  = h(x,y)  0 f(x,y ),  where  © denotes  convolution. 
An  ideal  device  would  have  h(x,y)  = 8(x,y).  The  Fourier  transform  /?,(£,  p)  of  the 
psf  is  generally  complex-valued  and  is  called  the  optical  transfer  function  (otf).  The 
absolute  value  of  the  otf  is  the  modulation  transfer  function  (mtf). 

In  [31],  it  is  noted  that  electron  optical  mtfs  are  often  nearly  Gaussian  in  shape, 
and  that  this  should  be  expected  from  the  Central  Limit  Theorem,  since  the  process 
of  converting  incoming  signal  photons  into  the  final  image  that  is  observed  on  a screen 
involves  many  intermediate  stages.  However,  it  is  also  observed  in  [31]  that  when  such 
mtfs  are  fitted  with  Gaussians,  the  fitted  curves  often  have  slimmer  tails  than  is  the 
case  in  the  true  mtfs. 

A systematic  study  of  electron  optical  mtf  measurements  is  the  subject  of  [17], 
[19],  and  [22].  There,  the  author  claims  the  empirical  discovery  that  a wide  variety  of 
electronic  imaging  devices,  including  phosphor  screens  and  some  types  of  photographic 
film,  have  otfs  that  are  well-described  by  (1)  with  1/2  < (3  < 1.  For  any 

given  device,  the  values  of  o-  and  [3  can  be  determined  using  specialized  graph  paper 
[23].  Other  instances  of  electron  optical  stable  laws  are  mentioned  in  [18],  [21],  and 
[25].  Analysis  of  the  physical  mechanisms  responsible  for  such  non  Gaussian  behavior 
is  not  included  in  these  works.  An  understanding  of  such  mechanisms  may  lead  to 
the  design  of  imaging  devices  with  low  values  of  [3.  The  latter  parameter  affects  the 
attenuation  of  high  frequency  information  in  the  recorded  image.  Deconvolution  of 
that  image  in  the  presence  of  noise,  is  generally  better  behaved  at  low  values  of  (3 
than  it  is  at  high  values  of  [3. 

The  characterization  (1)  is  useful  in  other  areas  of  optics.  The  otf  for  long- 
exposure  imaging  through  atmospheric  turbulence  [15],  is  known  to  be  given  by  (1) 
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with  (3  = 5/6,  and  a determined  by  atmospheric  conditions.  Also,  as  shown  in  [20], 
the  analytically  known  diffraction-limited  otf  for  a perfect  lens  [32,  p.  154],  can  be 
approximated  over  a wide  frequency  range  by  (1),  with  f3  = 3/4  and  a a properly 
chosen  function  of  the  cutoff  frequency. 

The  range  of  (3  values  discussed  above,  namely  1/2  < (3  < 1,  mirrors  that  found 
in  most  other  physical  contexts  where  symmetric  stable  laws  appear  or  are  surmised. 
Values  of  (3  <C  1/2  seem  to  be  relatively  rare  in  applications.  Examples  of  such  (3 
values  occur  in  [24],  where  mtf  data  for  56  different  kinds  of  photographic  film  are 
analyzed.  Good  agreement  is  found  when  these  data  are  fitted  with  (1),  and  the  pairs 
(n,  (3)  characterizing  each  of  these  56  mtfs  are  identified.  It  is  found  that  36  types  of 
film  have  mtfs  where  1/2  < (3  < 1.  The  remaining  20  types  have  mtfs  with  values  of 
j3  in  the  range  0.265  < f3  < 0.475. 

We  now  consider  imaging  systems  composed  of  various  elements  satisfying  (1). 
Such  systems  might  be  used  to  image  objects  through  a turbulent  atmosphere  or 
through  other  distorting  media  whose  oti’s  obey  (1).  The  resulting  composite  otf  has 
the  form 

(2)  /?,(£,  //)  = aA^+^)l3‘ , a,  >0,  0 < /3,  < 1. 

Such  an  object  corresponds  to  a multifractal  law  in  [4],  We  define  the  class  G to  be 
the  class  of  all  point  spread  functions  h(x,y)  with  Fourier  transforms  satisfying  (2). 
We  shall  be  interested  in  image  deblurring  problems 

(3)  Hf  = / h(x  - u,  y - v)f(u,  v)dudv  = h(x,  y ) <g>  f{x,  y)  = < fix,  y), 

Jr 2 

where  g(x,y)  is  the  recorded  blurred  image,  f(x,y ) is  the  desired  unblurred  image, 
and  h(x,y)  is  a known  point  spread  function  in  class  G.  The  blurred  image  g{x,y ) 
includes  noise,  which  is  viewed  as  a separate  additional  degradation, 

(-4)  fj(xcy)  = ge{x,y)  + n{x,y). 

Here,  ge(x/y)  is  the  blurred  image  that  would  have  been  recorded  in  the  absence  of 
noise,  and  n(x,y)  represents  the  cumulative  effects  of  all  errors  affecting  final  acqui- 
sition of  the  digitized  array  g(x,y).  The  unique  solution  of  (3)  when  the  right  hand 
side  is  ge(x,y),  is  the  exact  sharp  image  denoted  by  fe(x,y).  Thus 

(5)  h(x,y)  ® fe{x,y)  = ge{x,y). 

3.  Deblurring  with  the  SECB  method.  The  SECB  method  is  a direct  FFT- 
based  image  deblurring  technique  designed  for  equations  of  the  form  (3)  when  h(x,  y)  is 
known  and  satisfies  (2).  Theoretical  analysis  of  that  method,  along  with  error  bounds 
and  comparison  with  other  methods,  may  be  found  in  [5],  [6],  [7],  [8].  Significantly, 
the  method  does  not  impose  smoothness  constraints  on  the  unknown  image  f(x,y). 
For  class  G psfs,  we  may  define  fractional  powers  H1 , 0 < t < 1,  of  the  convolution 
integral  operator  H in  (3)  as  follows 

(6)  H'f  = {h'(G//)/(Gd)}  , 0 < t < 1. 

Class  G psfs  are  intimately  related  to  diffusion  processes,  and  solving  (3)  is  equivalent 
to  finding  the  initial  value  u(x,y,  0)  = f(x,y)  in  the  backwards  in  trine  problem  for 
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the  generalized  diffusion  equation 

U,  = - E/=1  A ,(-A)'-V  A,  = «,( 4tt2)-a,  I)  < f < 1. 

(7) 

"(i)  = u(x,y)- 

When  f{x,y)  is  known,  u(x,y,t)  = H1  f is  the  solution  of  (7)  at  time  t.  The  SECB 
method  is  a regularization  method  for  solving  the  ill-posed  problem  (7)  that  takes 
into  account  the  presence  of  noise  in  the  blurred  image  data  g(x,y)  at  / = 1.  The 
SECB  deblurred  image  /^  (,/•.  y)  is  obtained  in  closed  form  in  Fourier  space.  With  z 
denoting  the  complex  conjugate  of  c, 


(8) 


fH^n) 


M^v)g(Cv) 

|h(^7/)P+/v-->|l-/D'(e,7y)|2’ 


leading  to  p(x,y)  upon  inverse  transforming.  Here,  the  regularization  parameters 
A,  s are  positive  constants  that  depend  on  a-priori  information.  In  practice,  FFT 
algorithms  are  used  to  obtain  f'(x,y).  This  may  result  in  individual  pixel  values  that 
are  negative,  or  that  exceed  255,  the  maximum  value  in  an  8-bit  image.  Accordingly, 
all  negative  values  are  reset  to  the  value  zero,  and  all  values  exceeding  255  are  reset  to 
the  value  255.  For  512  x 512  images,  a single  trial  SECB  restoration  requires  about  5 
seconds  of  cpu  time  on  current  desktop  workstations.  We  may  also  form  and  display 

(9)  uHx,y,t)  = H'fHx,y), 


for  selected  decreasing  values  of  t lying  between  1 and  0.  This  simulates  marching 
backwards  in  time  in  (7),  and  allows  monitoring  the  gradual  deblurring  of  the  image. 
As  t f 0 the  partial  restorations  u^(x,y,t)  become  sharper.  However,  noise  and 
other  artifacts  typically  become  more  noticeable  as  t f 0.  Gradual  deblurring  allows 
detection  of  features  in  the  image  before  they  become  obscured  by  noise  or  ringing 
artifacts.  Such  marching  backwards  in  time  is  an  important  element  in  the  APEX 
method. 

It  should  be  noted  that  the  class  G is  only  a small  subclass  of  the  class  of  in- 
finitely divisible  densities,  [13].  The  latter  class  includes  multimodal  non  symmetric 
psfs,  associated  with  linear  diffusion  equations  more  complex  than  (7).  Detection  of 
such  psfs  from  blurred  image  data  would  require  considerable  extension  of  the  APEX 
method  discussed  below. 

4.  Non  uniqueness  in  blind  deconvolution.  Blind  deconvolution  of  images 
is  a mathematical  problem  that  is  not  fully  understood.  Well-documented  examples 
of  the  kinds  of  behavior  that  may  occur  are  of  particular  interest.  In  this  section, 
we  highlight  important  non  uniqueness  aspects  of  that  problem  that  are  helpful  in 
understanding  the  results  of  the  APEX  method.  Let  fe(x,y)  be  a given  exact  sharp 
image,  let  h(x,  y)  be  a Levy  point  spread  function,  and  let  g,.(x,  y)  — //.( x , y)<-  fe(x,  y). 
We  shall  show  that  given  the  blurred  image  gt.(x,y ),  there  are  in  general  many  point 
spread  functions  li,(x,y ) ^ h(x,y)  that  deblur  ge(x,y),  producing  high  quality  recon- 
structions fi{x,  y)  / /,  (./•,  //),  with  h ,(x.  y)  (■  f,(x.y)  « g,  (x,y). 

The  sharp  512  x 512  Sydney  image  /,  ( .r , y ) in  Figure  1(A)  was  synthetically 
blurred  by  convolution  with  a Cauchy  density  h(.r.y)  with  <*0  = 0.075,  fio  — 0.5. 
This  produced  the  blurred  image  g,  (.r.  y)  in  Figure  1(B).  To  avoid  distractions  caused 
by  noise,  the  blurred  image  ge(x,y)  in  this  experiment  was  computed  and  stored 
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Fig.  1.  Non  uniqueness  in  blind  deconvolution.  Distinct  point  spread  functions  exist  that 
produce  hiqh.  quality  reconstructions  from  same  blurred  image.  (A)  Original  sharp  512  x 512  Sydney 
image.  (B)  Synthetically  blurred  Sydney  image  created  by  convolution  with  Lorentzian  density  with 
qq  = 0.075,  /3o  = 0.5.  Blurred  image  computed  and  stored  in  64-bit  precision.  (C)  Deblurring 
of  image  (B)  using  correct  parameters  a = 0.075,  ft  = 0.5.  (D)  Deblurring  of  image  (B)  using 
a — 0.1301264,  ft  — 0.44298.  (E)  Deblurring  of  image  (B)  using  a = 0.1950345,  ft  = 0.403889.  (F) 
Deblurring  of  image  (B)  using  a — 0.2360994,  ft  = 0.369666.  Notice  that  images  (D),  (E),  and  (F) 
were  found  using  specific  pairs  (a,  ft)  where  a > ao  and  ft  < fto-  All  deblurred  images  obtained 
using  SECB  procedure  with  s = 0.001  and  h — 10000. 
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iii  64-bit,  precision.  Deblurring  this  noiseless  image  with  the  correct  psf  values  a = 
0.075,  /3  = 0.5,  produces  Figure  1(C).  This  is  in  excellent  visual  agreement  with 
fe(x,y)  m Figure  1(A),  as  expected.  However,  the  visual  quality  in  Figures  1(D), 
1(E).  and  1(F)  is  generally  as  good  as  that  in  Figure  1(C).  The  latter  three  images 
were  deblurred  with  Levy  densities  with  values  (q-,/1)  where  o > ci0,  /j  < /io , and 
they  differ  from  Figure  1(A)  in  contrast  and  brightness.  All  deblurred  images  were 
obtained  using  the  SECB  method  with  s = 0.001  and  I\  = 10000.  One  dimensional 
cross  sections  of  the  four  distinct  psfs  used  in  Figure  1 are  displayed  in  Figure  2. 
These  psfs  also  exhibit  distinct  heavy  tail  behavior  not  shown  in  Figure  2. 

One  can  imagine  four  photographers,  simultaneously  photographing  the  identical 
scene  depicted  in  Figure  1(A),  yet  producing  the  four  distinct  images  shown  in  Figures 
1(C),  1(D),  1(E),  and  1(F),  through  use  of  different  lenses,  film,  filters,  exposures, 
printing,  and  the  like.  In  practice,  given  only  the  blurred  image  in  Figure'  1(B),  any 
one  of  these  four  restorations  would  be  considered  highly  successful.  Convolution  of 
each  reconstruction  with  its  corresponding  psf  in  Figure  2.  reproduces  the  blurred 
image  in  Figure  1(B). 

For  any  restoration  /(./•.  //)  of  the  exact  image  fe(x,y)  in  Figure  1(A),  anti  any 
norm  ||  ||,  we  can  evaluate  the  relative  error  ||  / — /,  ||  / ||  fe  ||.  Define  the  discrete 

Id,  L J . and  Hm  norms  as  follows 

(10)  11/11-2=  {a'-2EC=.I/(-'2  9)P}'/', 

ii  / ii*  - = { iv-2  e£k,=.  ( i + e + tY"  \m.  -i)P  }1/2 . 

The  relative  errors  in  the  L] , L1 . II 1 and  //  ’ norms,  for  each  of  the  four  restora- 
tions in  Figure  1,  are  shown  in  Table  1.  As  might  be  expected,  image  (C)  is  the  closest 
to  image  (A)  in  each  of  these  norms,  since  the  correct  psf  values  were  used  to  obtain 
image  (C)  from  image  (B).  It  is  also  evident  from  Table  1 that  the  four  restorations 
are  distinct  from  one  another,  since  they  differ  from  image  (A)  by  different  amounts. 
Most  important,  the  fact  that  image  (E)  is  a significantly  poorer  approximation  to 
image  (A)  in  these  norms  than  is  image  (C),  does  not  imply  that  image  (E)  is  an 
inaccurate  representation  of  the  visual  scene  depicted  in  image  (A).  Notice  also  that 
image  (F)  is  not  as  sharp  as  image  (E),  although  it  is  closer  to  image  (A)  in  three  of 
the  four  norms. 

Iterative  algorithms  are  the  most  common  approach  to  blind  deconvolution.  Con- 
vergence proofs  for  such  iterative  procedures  are  seldom  available.  The  above  example 
illustrates  some  of  the  difficulties  underlying  any  analysis  of  convergence.  Such  anal- 
ysis should  allow  for  the  possibility  of  infinitely  many  useful  limit  points,  while  the 
mathematical  characterization  of  such  limit  points  is  not  obvious.  Moreover,  as  is 
evident  from  Table  1 and  has  been  known  for  some  time,  the  use  of  Lp  or  H‘"  norms 
in  assessing  the  visual  quality  of  a reconstruction  can  be  misleading. 

5.  Marching  backwards  in  time  and  the  APEX  method.  The  APEX 
method  is  a blind  deconvolution  technique  based  on  detecting  class  G psf  signat  ures 
by  appropriate  1-D  Fourier  analysis  of  the  blurred  image  y(.i\  y).  The  detected  psf 
parameters  are  then  input  into  the  SECB  algorithm  to  deblur  the  image.  Let  /,.(.;■.  //) 
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FOUR  DISTINCT  PSFS  THAT  DEBLUR  SYDNEY  IMAGE 


SPATIAL  VARIABLE  X (IN  PIXELS) 


C5 


Os 


Fig.  2.  Four  distinct  point  spread  functions  that  deblur  image  (B)  in  Figure  1.  Curves  C , D,  E 
and  F are  l-D  cross  sections  of  the  512  x 512  psfs  that  respectively  produced  images  (C),  (D),  (E) 
and  (F)  in  Figure  1.  These  psfs  also  exhibit  distinct  heavy  tail  behavior. 

TABLE  1 


Relative  errors  in  various  norms  for  the  four  deblurred  images  in  Figure  1. 


Restoration 

Parameters  o,  (3 

Ll 

L 2 

Hl 

H 5 

Image  (C) 

a = 0.075,  /3  = 0.500 

2.13  % 

3.52  % 

4.13  % 

19.66  % 

Image  (D) 

a = 0.130,  f3  = 0.443 

6.63  % 

8.37  % 

8.67  % 

21.11  % 

Image  (E) 

a = 0.195,  /3  = 0.404 

12.64  % 

15.53  % 

15.75  % 

25.52  % 

Image  (F) 

a = 0.236,  f3  — 0.370 

12.54  % 

15.08  % 

15.31  % 

26.17  % 

be  an  exact  sharp  image  as  in  (5).  Since  fe(x,y)  > 0 

(11)  |/e (f , v)\  < / fe (x,  y)dxdy  = fe (0, 0)  = a > 0. 

Also,  since  ge(x,y)  = h.(x,y ) 0 fe(x,y)  and  h(x,y)  is  a probability  density, 

(12)  ge (0, 0)  — / ge(xpy)dxdy  = I f,  (x.y)dxdg  = /,  ((), 0)  = rr  > 0. 

Jr 2 Jr 2 

Using  a as  a normalizing  constant,  we  may  normalize  Fourier  transform  quantities 
q{ev)  dividing  by  a.  Let 


(13) 
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Fig.  3.  Behavior  of  normalized  Fourier  transform  in  types  of  blurred  images  g(x,y)  considered 
in  present  paper.  (A)  log  \g*  (£,  0)|  in  Cindy  Crawford  image.  (B)  log|r)*(£.0)|  in  Washington  DC 
Landsat  image.  While  local  behavior  is  highly  oscillatory,  global  behavior  is  generally  monotone 
decreasing  and  convex. 


denote  the  normalized  quantity.  The  function  |/e  (£, rf)\  is  highly  oscillatory,  and 
- * 

0 < l/e  | < I-  Since  fe(x,y)  is  real,  its  Fourier  transform  is  conjugate  symmetric. 

» * 

Therefore,  the  function  \ fe  (£,'//)  | is  symmetric  about,  the  origin  on  any  line  through 
the  origin  in  the  (£,77)  plane.  The  same  is  true  for  the  blurred  image  data  \g*(£,  7/ ) | . 

All  blurred  images  in  this  and  the  next  section  are  of  size  512  x 512  and  quantized 
at  8-bits  per  pixel.  For  any  blurred  image  g(x,y),  the  discrete  Fourier  transform  is  a 
512  x 512  array  of  complex  numbers,  which  we  again  denote  by  </(£,  77)  for  simplicity. 
The  ‘frequencies’  £,?/  are  now  integers  lying  between  —256  and  256,  and  the  zero 
frequency  is  at  the  center  of  the  transform  array.  This  ordering  is  achieved  by  pre- 
multiplying g(x,  y)  by  ( — l)J  + y.  We  shall  be  interested  in  the  values  of  such  transforms 
along  single  lines  through  the  origin.  The  discrete  transforms  |<7*(£,  0)|,  and  |<y*(0, 77) | 
are  immediately  available.  Image  rotation  may  be  used  to  obtain  discrete  transforms 
along  other  directions.  All  1-D  Fourier  domain  plots  shown  in  this  paper  are  taken 
along  the  axis  7/  = 0 in  the  (£,77)  plane.  In  these  plots,  the  zero  frequency  is  at  the 
center  of  the  horizontal  axis,  and  the  graphs  are  necessarily  symmetric  about  the 
vertical  line  £ = 0.  Examples  of  such  plots  are  shown  in  Figures  3,  5,  and  10. 

The  class  of  blurred  images  g{x,  y)  considered  in  the  present  paper  can  be  de- 
scribed in  terms  of  the  behavior  of  log  \g*(£,T))\  along  lines  through  the  origin  in  the 
(£,77)  plane.  While  local  behavior  is  highly  oscillatory,  global  behavior  is  generally 
monotone  decreasing  and  convex.  This  is  shown  in  Figure  3 for  two  typical  images, 
along  the  line  77  = 0.  I11  [9],  a large  class  of  images  with  that  property  was  exhibited, 
the  class  W.  The  blurred  images  considered  here  may  be  loosely  characterized  as 
being  in  class  W.  Not  all  blurred  images  may  be  so  characterized.  For  example,  if 
the  Cindy  Crawford  image  g(x,y)  in  Figure  3(A)  were  convolved  with  a,  wide  Gaus- 
sian psf  to  form  a new  blurred  image  gi(x,  y),  global  behavior  in  log  |</i*(£,0)|, 
away  from  the  origin,  would  be  monotone  decreasing  and  concave.  Application  of  the 
APEX  method  to  several  concave  examples  is  discussed  in  [9].  Convolution  of  Figure 
3(A)  with  a.  defocus  psf  produces  a different  kind  of  blurred  image  //■_> (;/\  ;</),  and  global 
behavior  in  log  |fyV(£,  0)|  is  neither  concave  nor  convex.  Instead,  there  is  a regular 
pattern  of  sharp  singularities  corresponding  to  successive  zeroes  of  the  defocus  off. 
Use  of  the  APEX  method  in  the  manner  to  be  described  below,  is  intended  only  for 
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blurred  images  with  Fourier  behavior  analogous  to  that  shown  in  Figure  3. 

The  APEX  method  is  based  on  the  following  observations.  In  the  basic  relation 

(14)  g(*,y)  = h{x,y)  ® fe(x,y)  + n(x,y), 
we  may  safely  assume  that  the  noise  ii(x,y)  satisfies 

(15)  / \n(x,y)\dxdy  / f,.(x,y)dxdy  = a > 0, 

Jr2  ■)  r- 

so  that, 

(16)  |b*(£,r/)|  <C  1. 

Consider  the  case  where  the  otf  is  a pure  Levy  density  />(£,?/)  = e~a^  +n  1 . Since 
!)  = 9e  + II 

(17)  log|<nC'/y)|  =\og\e-a{?  + ,'2),i  l*(Z,v)  +n'(Z,ii)\. 

Let  = {(£,?/)  | £2  + T)2  < ur}  be  a neighborhood  of  the  origin  where 

(18)  e>-a(€WC|//(^7/)|»|7r(4i7?)|. 

Such  an  il  exists  since  (18)  is  true  for  £ = ?/  = 0 in  view  of  (16).  The  radius  u > 0 of 
11  decreases  as  a and  n increase.  For  (£,//)  £ 11  we  have 

(19)  log|<T(£,t?)|  » -o(£2  + rfY  + log  |/e*(£, r;)|. 

Because  of  the  radial  symmetry  in  the  psf,  it  is  sufficient  to  consider  (19)  along  a 
single  line  through  the  origin  in  the  (£,??)  plane.  Choosing  the  line  y = 0,  we  have 

(20)  log  l^*(C  0)|  ~ — o|£|2/i  + log  |/e*(£,  0)|,  kc|<u,’. 

Some  type  of  a-priori  information  about  fe(x,y)  is  necessary  for  blind  deconvolu- 
tion. In  (20),  knowledge  of  log|/e  (£,0)1  on  |£|  < uj  would  immediately  yield  Q'|£|2'^ 
on  that  interval.  Moreover,  any  other  line  through  the  origin  could  have  been  used  in 
(19).  However,  such  detailed  knowledge  is  unlikely  in  practice.  The  APEX  method 
seeks  to  identify  a useful  psf  from  (20)  without  prior  knowledge  of  log  |/e  (£,  0)|.  The 
method  assumes  instead  that  fe(x,  y)  is  a recognizable  object,  and  typically  requires 
several  interactive  trials  before  locating  a suitable  psf.  As  previously  noted,  such 
trial  SECB  restorations  are  easily  obtained.  Here,  prior  information  about  fe(x,y) 
is  disguised  in  the  form  of  user  recognition  or  rejection  of  the  restored  image,  and 
that  constraint  is  applied  at  the  end  of  the  reconstruction  phase,  rather  than  at  the 
beginning  of  the  detection  phase. 

~ * 

In  the  absence  of  information  about  log|/e  (£,0)|,  we  replace  it  by  a negative 
constant  —.4  in  (20).  For  any  A > 0,  the  approximation 

(21)  log|r(C0)|  w -a\Z\2l3-A, 

is  not  valid  near  £ = 0,  since  the  curve  u(£)  = — o|£|2 ,H  — A,  has  —A  as  its  apex. 
Choosing  a value  of  .4  > 0,  we  best  fit  log  |<y*(£,0)|  with  u(£)  — — oj£|2tf  — .4  on  the 
interval  |£|  < a;,  using  nonlinear  least  squares  algorithms.  The  resulting  fit  is  close 
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Fig.  4.  Enhancement  of  Cindy  Crawford  image  by  marching  backwards  from  t = 1 with 
APEX  detected  psf.  Image  sequence  shows  gradual  increase  in  contrast  as  I decreases.  Undesirable 
artifacts  at  t = 0 indicate  continuation  backwards  in  time  has  proceeded  too  far.  Best  results  are 
highly  subjective  m this  case,  but  probably  occur  at  some  I > 0.5.  Note  sharpness  of  earrings  near 
t = 0.5. 


only  for  £ away  from  the  origin.  The  returned  values  for  a and  ft  are  then  used  in 
the  SECB  deblurring  algorithm.  Different  values  of  .4  return  different  pairs  (cv,  (3). 
Experience  indicates  that  useful  values  of  .4  generally  lie  in  the  interval  2 < A < 6. 
Increasing  the  value  of  A decreases  the  curvature  of  u{£)  at  £ = 0,  resulting  in  a 
larger  value  of  ft  together  with  a smaller  value  of  a.  A value  of  .4  > 0 that  returns 
ft  > 1 is  clearly  too  large,  as  ft  > I is  impossible  for  probability  density  functions 
[13].  Decreasing  A has  the  opposite  effect,  producing  lower  values  of  ft  and  higher 
values  of  a.  As  a rule,  deconvolution  is  better  behaved  at  lower  values  of  ft  than  it  is 
when  ft  « 1.  A significant  observation  is  that  an  linage  blurred  with  a pair  (ao,fto) 
can  often  be  successfully  deblurred  with  an  appropriate  pair  (a,  ft),  where  o > «o 
and  ft  < ft0  . Examples  of  this  phenomenon  were  shown  in  Figure  1 in  connection 
with  the  blurred  Sydney  image.  An  effective  interactive  framework  for  performing  the 
above  least  squares  fitting  is  the  fit  command  in  DATAPLOT  [16].  This  is  a high- 
level  English-syntax  graphics  arid  analysis  software  package  developed  at  the  National 
Institute  of  Standards  and  Technology.  This  software  tool  was  used  throughout  this 
paper. 

The  following  version  of  the  APEX  method,  using  the  SECB  inarching  backwards 
in  time  option  (9),  has  been  found  useful  in  a variety  of  image  enhancement  problems 
where  the  image  g(x,  y)  is  such  that  log  0)|  is  generally  globally  monotone 
decreasing  and  convex,  as  shown  in  Figure  3.  Choose  a value  of  .4  > 0 in  (21)  such 
that  the  least  squares  fit  develops  a slight  cusp  at  f = 0.  Using  the  returned  pair 
(a,  ft)  in  the  SECB  method,  obtain  a.  sequence  u'(x,y,t)  of  partial  restorations  as 
t decreases  from  t = 1,  as  illustrated  in  the  Cindy  Crawford  sequence  in  Figure  4. 
With  a good  choice  of  A,  high  quality  restorations  will  be  found  at  positive  values  of 
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t,  and  these  will  gradually  deteriorate  as  t 4 0.  Typically,  the  restoration  at  t = 0 
will  exhibit  undesirable  artifacts,  indicating  that  continuation  backwards  in  time  has 
proceeded  too  far  in  (7).  Terminating  the  continuation  at  some  appropriate  t = t.\  > 0, 
is  equivalent  to  rescaling  the  value  of  a without  changing  the  value  of  (3.  If  the  pair 
(a,/i)  produces  a high  quality  restoration  at  t.  = t\  > 0,  the  pair  (aq,/3),  where 
aj  = (1  - t i)cv,  will  produce  the  same  quality  results  at  t = 0.  In  general,  there  will 
be  many  values  of  A in  (21)  returning  pairs  {a,  (3)  that  produce  good  reconstructions 
at  some  tap  > 0.  A large  number  of  distinct  pairs  («*,/4*)  can  thus  be  found  that 
produce  useful,  but  distinct,  results  at  t = 0.  Indeed,  this  is  the  process  that  was  used 
to  obtain  the  four  psfs  shown  in  Figure  2. 

We  have  been  assuming  //.(£,//)  to  be  a pure  Levy  otf  in  (14).  For  more  general 
class  G otfs  (2),  we  may  still  use  the  approximation  log|<?*(£, 0)|  ~ — a|£|2/^  - .4, 
and  apply  the  same  technique  to  extract  a suitable  pair  (cv,/l)  from  the  blurred  image. 
Here,  the  returned  APEX  values  may  be  considered  average  values  for  the  o,,  (3,  in 
(2).  producing  a single  pure  Levy  otf  approximating  the  composite  otf. 

6.  Application  to  real  images.  The  developments  in  sections  2 through  5 
are  predicated  on  two  assumptions.  The  first  assumption  is  that  the  blurred  image 
g(.r,  y)  obeys  the  simple  convolution  equation  (3)  rather  than  a more  general,  possibly 
nonlinear,  integral  equation 


In  addition  to  linearity,  (3)  implies  that,  the  blur  is  isoplanatic.  The  second  assumption 
is  that  the  point  spread  function  h(x,y)  belongs  to  a restricted  class  of  unimodai, 
radially  symmetric,  probability  density  functions,  the  class  G defined  in  (2).  In  [9], 
successful  blind  deconvolution  of  synthetically  blurred  images,  with  added  noise,  was 
demonstrated.  Such  synthetically  blurred  images  necessarily  obey  (2)  and  (3). 

The  applicability  of  the  preceding  theory  to  real  blurred  images  is  by  no  means 
assured.  Deviations  from  linearity,  isoplanatism,  unimodality,  and  radial  symmetry 
are  possible.  Moreover,  the  class  G excludes  motion  and  defocus  blurs.  In  addition, 
the  types  and  intensities  of  noise  processes  in  real  images  may  differ  fundamentally 
from  the  noise  models  typically  used  in  numerical  experiments.  Therefore,  only  limited 
success  on  a narrow  class  of  images  can  be  expected  in  real  applications. 

The  examples  discussed  below  involve  images  obtained  from  multiple  sources  using 
diverse  imaging  modalities.  Some  of  these  images  have  been  used  as  test  images  in 
the  literature.  In  this  paper,  each  of  these  images  is  assumed  to  have  been  blurred  by 
some  unknown  process,  and  we  seek  to  improve  visual  quality  by  APEX  processing. 
All  images  are  of  size  512  x 512  and  quantized  at  8 bits  per  pixel. 

Our  first,  example  is  a well-known  English  village  image  denoted  by  g(x,y),  and 
shown  in  Figure  5(A)  together  with  log  |ry*(£,  0)|  on  |£|  < 250.  The  plot,  displays 
globally  convex  monotone  behavior.  In  Figure  5(B),  the  APEX  fit  of  log|g*(£, 0)| 
with  u(f)  = — o|£|2/3  — .4,  on  the  interval  |£|  < 200,  is  shown.  With  .4  = 3.75,  the  fit 
develops  a cusp  at  £ = 0 and  returns  o = 0.251274,  (3  = 0.242246.  With  these  psf 
parameters,  SECB  deblurring  using  s = 0.01,  A = 1300,  and  continuation  backwards 
in  time  terminated  at.  t — 0.5,  produces  Figure  6(B).  This  is  compared  with  the 
original  in  Figure  6(A). 

The  extent  of  sharpening  in  Figure  6(B)  becomes  evident  when  zooming  on  se- 
lected parts  of  the  image.  In  Figure  7,  rooflines  on  the  first  three  houses  are  compared 
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Fig.  5.  APEX  method  of  psf  detection.  (A)  log  |g*  (£,  0)|  on  |£|  < 250  m 8-bit  English,  village 
image.  (B)  Least  squares  fit  of  log  | g*  (f , 0)|  with  ;/,(£)  = —a  |£|2^  —3.75  on  |£|  < 200.  develops  cusp 
at  £ = 0 and  returns  a = 0.251274,  /3  — 0.242246. 


A 


B 


Fig.  6.  Enhancement  of  English  village  image.  (A)  Original  8-bit  image.  (B)  SECB  deblurred 
image  using  s — 0.01,  l\  = 1300,  with  APEX  detected  values  a = 0.251274,  /l  = 0.242246.  and 
continuation  backwards  in  time  terminated  at  t = 0.5. 


before  and  after  APEX  processing.  There  is  noticeable  enhancement  of  structural  de- 
tail in  the  roof  shingles  and  stone  fronts  of  the  three  houses  in  Figure  7(B).  In  Figure 
8(B),  Holstein  cows  grazing  in  the  meadow,  not  previously  identifiable,  are  clearly 
visible.  So  are  individual  chimney  bricks.  In  Figure  9(B),  buildings  in  the  distance, 
not  readily  noticed  in  Figure  9(A),  become  well-defined. 

It  should  be  noted  that  use  of  a different  value  of  .4,  and/or  a different  neigh- 
borhood of  the  origin  H in  Figure  5(B),  may  return  a different  psf  pair  (cv,/J).  In 
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Fig.  8.  Extent  of  sharpening  in  English  village  scene  becomes  more  evident  with  zooming. 
Holstein  cows  grazing  in  m.eadow  in  image  (B)  are  not  readily  identifiable  in  image  (A). 


that  case,  backwards  continuation  in  the  SECB  method  may  need  to  be  terminated 
at  some  other  value  of  t to  obtain  the  best  image.  However,  with  good  choices  of  A 
and  {l,  the  latter  image  would  again  be  a high  quality  representation  of  the  visual 
scene  in  Figure  6(B),  while  differing  from  Figure  6(B)  at  individual  pixels.  This  is 
the  non  uniqueness  phenomenon  previously  discussed  in  connection  with  the  Sydney 
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Fig.  9.  Extent  of  sharpening  m English,  village  image  becomes  more  evident  with  zooming. 
Enhanced  image  (B)  shows  buildings  in  distance  not  immediately  apparent  in  original  image  (A). 
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FOURIER  TRANSFORM  BEHAVIOR  IS  ENGLISH  VILLAGE  IMAGE 
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Fig.  10.  APEX  processing  significantly  alters  Fourier  transform  behavior.  (A)  English  village 
image  before  and  after  processing.  (B)  F15  terrain  image  in  Figure  12  before  and  after  processing. 
Behavior  shown  in  (B)  is  exceptional.  All  other  examples  in  present  paper  conform  with,  behavior 
shown  in  (A). 


image  in  Figure  1. 

Deconvolution  of  Figure  6(A)  with  the  above  APEX-detected  psf  significantly 
alters  its  Fourier  transform.  As  shown  in  Figure  10(A),  the  Fourier  transform  in 
Figure  6(B)  (dashed  curve),  decays  less  rapidly  as  |£|  increases  than  was  the  case  in 
the  original  Figure  6(A)  (solid  curve).  Evidently,  APEX  processing  amplifies  high 
frequency  image  components  in  a stable  coherent  fashion,  resulting  in  the  overall 
improvements  visible  in  Figures  6 through  9.  The  ‘before  and  after’  Fourier  transform 
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Fig.  11.  Enhancement  of  boat  image.  APEX  method  with  A = 4.0  on  |£|  < 250,  yields 
a = 0.518155,  = 0.215083.  Using  these  parameters  with  s = 0.01,  K = 1300,  and  backwards 

continuation  terminated  at  t = 0.5,  SECB  method  applied  to  image  (A)  produces  image  (B).  Number 
7 2 7 on  side  of  boat  in  image  (B)  was  not  easily  identifiable  in  image  (A). 


Fig.  12.  Striking  enhancement  of  terrain  features  in  FI5  image.  APEX  method  with  A = 3.5 
on  |£|  < 250,  yields  a = 0.85G096,  ft  = 0.107289.  Using  these  parameters  with  s = 0.01,  A = 1000, 
and  backwards  continuation  terminated  at  t — 0.25,  SECB  method  applied  to  image  (A)  produced 
image  (B).  Condensation  trails  behind  aircraft  in  image  (B)  not  immediately  evident  m image  (A). 
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Fig.  13.  Enhancement  of  Washington  DC  Landsat  image.  APEX  method  with  A = 4.25  on 
|£|  < 250,  yields  a = 0.540825,  f3  = 0.182410.  Using  these  parameters  with  s = 0.01,  K = 1300,  and 
backwards  continuation  terminated  at  t = 0.5,  SECB  method  applied  to  image  (A)  produced  image 
(B).  Increased  resolution  in  image  (B)  improves  definition  of  several  landmarks  and  thoroughfares. 


A B 


Fig.  14.  Enhancement  of  scanning  electron  microscope  image  of  a mosquito’s  head  showing 
compound  eye.  APEX  method  with  A = 4.0  on  |£|  < 250,  yields  a = 0.734259,  (3  = 0.156963.  Using 
these  parameters  with  s = 0.001,  K = 10,  and  backwards  continuation  terminated  at  t = 0.4,  SECB 
method  applied  to  image  (A)  produced  image  (B).  APEX  processing  enhances  contrast  and  brings 
eye  into  sharper  focus. 
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Fig.  15.  Enhancement  of  sagittal  MRI  brain  image.  APEX  method  with  A — 4.0  on  |£|  < 250, 
yields  a = 0.333267,  (3  — 0.209416.  Using  these  parameters  with  s = 0.01,  K = 1300,  and  backwards 
continuation  terminated  at  t = 0.35,  SECB  procedure  applied  to  image  (A)  produced  image  (B). 
APEX  processing  noticeably  improves  feature  definition  in  areas  between  two  and  four  o’clock. 


Fig.  16.  Enhancement  of  transverse  PET  brain  image.  APEX  method  with  A = 5.0  on 
|£|  < 250,  yields  a = 0.198931,  (3  = 0.284449.  Using  these  parameters  with  s = 0.001,  K = 5.0, 
and  backwards  continuation  terminated  at  t = 0.6,  SECB  procedure  applied  to  image  (A)  produced 
image  (B).  Bright  spots  in  enhanced  image  (B),  indicating  areas  of  the  brain  responding  to  applied 
external  stimuli,  are  barely  visible  m original  image  (A). 
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Fig.  17.  Enhancement  of  Whirlpool  Galaxy  (M51)  image.  APEX  method  with  A = 4.0  on 
|£|  < 250,  yields  a = 0.451615,  f3  = 0.221955.  Using  these  parameters  with  s = 0.001,  A’  = 5.0. 
and  backwards  continuation  terminated  at  t = 0.5,  SECB  applied  to  image  (A)  produced  image  (B). 
APEX  processing  increases  resolution  and  enhances  luminosity  in  spiral  arms  and  galactic  cores. 


pattern  shown  in  Figure  10(A)  occurs  in  every  example  discussed  in  this  paper,  with 
the  exception  of  the  F15  image  in  Figure  12.  The  anomalous  behavior  in  that  case  is 
shown  in  Figure  10(B). 

The  next  example  is  the  boat  image  in  Figure  11(A).  With  A — 4.0,  the  APEX  fit 
on  |£|  < 250  returned  o = 0.518155,  (3  = 0.215083.  Using  these  values  in  the  SECB 
method,  with  s = 0.01,  A = 1300,  and  continuation  terminated  at  t = 0.5,  produced 
Figure  11(B).  Enhancement  has  now  rendered  visible  the  number  7 2 7 on  the  left 
side  of  the  boat.  Other  identifiable  details  include  the  stripe  along  the  left  trouser 
leg  of  the  man  on  the  ground,  the  lettering  on  the  sign  hanging  from  the  boat,  to  his 
right,  and  part  of  the  stone  work  and  stairway  to  the  left  of  the  lighthouse. 

The  F15  plane  image  in  Figure  12(A)  is  another  interesting  example.  The  aim 
here  is  to  enhance  the  background  terrain.  With  A = 3.5,  the  APEX  fit  on  |^|  < 250 
develops  a cusp  at  £ = 0 and  returns  a = 0.856096,  (3  = 0.107289.  Using  these 
values  in  the  SECB  method,  with  s = 0.01,  I\  = 1000,  and  backwards  continuation 
terminated  at  t = 0.25,  produces  rather  striking  enhancement  of  the  ground  features 
in  Figure  12(B).  This  example  is  noteworthy  on  two  counts:  the  exceptionally  low 
value  of  f3  detected  by  the  APEX  method,  and  the  previously  mentioned  unexpected 
Fourier  behavior  shown  in  Figure  10(B). 

Beginning  with  Figure  1,  all  of  the  examples  discussed  so  far  involve  images  of 
familiar  objects.  This  allows  for  relatively  easy  evaluation  of  the  results  of  APEX 
processing.  The  next  five  examples  involve  less  familiar  objects.  Moreover,  fine  de- 
tails visible  on  a modern  high  resolution  computer  screen  are  sometimes  lost  in  t he 
printing  process.  Consequently,  improvements  in  image  quality  in  some  of  the  next 
examples  may  seem  less  obvious  than  in  previous  examples.  At  the  same  time,  the 
performance  of  the  APEX  method  in  reconstructing  real  details  of  familiar  objects, 
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provides  a.  measure  of  confidence  in  the  results  obtained  when  that  method  is  applied 
to  unfamiliar  objects. 

Figure  13(A)  is  a Landsat  image  of  the  Washington  DC  area.  With  A = 4.25,  the 
APEX  fit  on  |f  | < 250  returns  a = 0.540825,  f3  = 0.182410.  Using  these  parameters  in 
the  SECB  method  with  s = 0.01,  K — 1300,  and  continuation  terminated  at  t = 0.5, 
produces  Figure  13(B).  There  is  a significant  increase  in  resolution  in  Figure  13(B), 
which  improves  definition  of  several  landmarks  and  thoroughfares.  The  Washington 
Monument,  the  bridges  over  the  Potomac,  Pennsylvannia  and  Maryland  Avenues 
radiating  from  the  Capitol,  Massachusetts  Avenue  to  the  north,  and  Virginia  Avenue 
and  the  Southeast  Freeway  to  the  south,  are  some  of  the  features  that  are  more  easily 
identified  in  the  enhanced  image. 

Figure  14(A)  is  a scanning  electron  microscope  image  of  a mosquito’s  head.  A 
prominent  feature  is  the  insect’s  compound  eye.  With  .4  = 4.0,  the  APEX  fit  on  |f  | < 
250  yields  a = 0.734259,  f3  = 0.15G9G3.  Using  these  values  in  the  SECB  method  with 
s = 0.001,  I\  — 10.0,  and  backwards  continuation  terminated  at  t — 0.4,  produces 
Figure  14(B).  Evidently,  APEX  processing  results  in  significant  overall  improvement. 
In  particular,  the  eye  appears  in  much  sharper  focus. 

The  sagittal  MRI  brain  image  in  Figure  15(A)  has  been  used  as  a test  sharp  image 
in  previous  publications.  In  [5]  and  [7],  synthetically  blurred  versions  of  that  sharp 
image  were  used  in  a comparative  evaluation  of  restoration  algorithms  when  the  psf  is 
known.  Here,  we  consider  further  sharpening  the  sharp  image  by  blind  deconvolution. 
With  .4  = 4.0.  the  APEX  fit  on  |f|  < 250  returns  a — 0.333267,  [3  = 0.209416.  Using 
these  parameters  in  the  SECB  procedure  with  s = 0.01,  K = 1300,  and  continuation 
terminated  at,  t = 0.35,  produced  the  image  in  Figure  15(B).  Substantial  improvement 
is  apparent  over  the  whole  image.  In  the  sector  between  two  and  four  o’clock  in 
particular,  sharpening  of  structural  detail  significantly  improves  feature  definition. 

In  PET  imaging,  a positron  emitting  radionuclide  is  injected  into  the  patient  and 
used  to  tag  glucose  molecules  in  their  course  through  the  brain.  The  metabolic  rate 
of  glucose  is  a key  parameter  that  measures  cerebral  function  and  reflects  the  extent 
to  which  regions  of  the  brain  are  working  or  failing  to  work.  Performing  specific 
mental  tasks  activates  various  parts  of  the  brain,  causing  increased  glucose  uptake 
and  hence  increased  positron  emission.  Centers  of  activity  translate  into  relatively 
bright  spots  in  the  PET  image.  However,  blurring  by  the  scanner  psf  tends  to  average 
out  such  relative  differences,  resulting  in  loss  of  contrast.  Figure  16(A)  is  a PET 
image  of  a transverse  slice  through  the  brain.  Blind  deconvolution  is  used  to  enhance 
that  image.  With  .4  = 5.0,  the  APEX  fit  on  |£|  < 250  returns  cv  = 0.198931,  (3  = 
0.284449.  Using  these  parameters  in  the  SECB  method  with  s = 0.001,  I\  = 5.0, 
and  backwards  continuation  terminated  at  t = 0.6.  produces  Figure  16(B).  Note  that 
both  images  in  Figure  16  show  identical  features,  but  contrast  has  been  increased  in 
the  APEX  processed  image,  with  some  regions  becoming  darker  while  others  have 
become  lighter.  In  particular,  several  bright  spots  appear  in  Figure  16(B)  that  were 
not  readily  apparent  in  the  original  image. 

Our  last  example  is  the  Whirlpool  galaxy  (M51)  in  Figure  17(A).  With  .4  = 4.0, 
the  APEX  fit  on  |£|  < 250  yields  a = 0.451615,  /3  = 0.221955.  Using  these  values  in 
the  SECB  method  with  s = 0.001,  I\  = 5.0,  and  backwards  continuation  terminated 
at  t = 0.5,  produced  Figure  17(B).  In  the  enhanced  image,  the  spiral  arms  are  more 
luminous  and  better  defined,  and  the  luminous  cores  are  larger  in  both  the  spiral 
galaxy  and  its  companion.  The  dark  connecting  arm  between  the  two  galaxies  is  also 
more  clearly  defined.  These  enhancements  are  due  to  a change  in  Fourier  transform 
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behavior  brought  about  by  deconvolution  with  the  APEX-detected  psf.  This  change 
in  Fourier  behavior  is  similar  to  that  shown  in  Figure  10(A).  although  it  is  more 
pronounced.  A concomitant  effect  of  deconvolution  is  amplification  of  data  noise, 
which  now  becomes  visible  against  the  dark  background  in  Figure  17(B). 

Clearly,  in  this  galaxy  image  as  in  the  preceding  PET  image,  there  is  no  way 
of  knowing  whether  or  not  the  enhanced  image  conforms  with  reality.  Conceivably, 
the  increased  luminosity  in  Figure  17(B)  may  be  exaggerated.  However,  the  bright 
areas  along  the  galactic  arms  in  Figure  17(B),  as  well  as  the  bright  spots  in  Figure 
1G(B),  did  not  materialize  spontaneously.  These  areas  must  have  been  just  below 
some  brightness  threshold  in  the  original  image,  and  APEX  processing  has  served  the 
very  useful  purpose  of  revealing  their  presence.  If  such  areas  appear  overenhanced, 
this  can  be  corrected  by  repeating  the  SECB  procedure  and  terminating  continuation 
at  higher  values  of  / . 

7.  Concluding  remarks.  Setting  aside  all  theoretical  considerations,  a practi- 
cal enhancement  technique  has  been  presented  that  can  sharpen  significant  classes  of 
images,  originating  from  diverse  imaging  modalities.  One  important  feature  of  the 
above  approach  is  its  fast  implementation  on  desktop  platforms.  Even  with  large  size 
images,  numerous  trial  restorations  can  be  accomplished  in  a few  minutes  of  cpu  time. 
This  makes  for  easy  hue  tuning  of  parameters.  Whether  or  not  APEX  processing  sig- 
nificantly improves  a given  image  can  generally  be  quickly  decided.  Once  improvement 
is  detected,  fine  tuning  must  be  used  to  obtain  optimal  results.  Here,  another  impor- 
tant feature  of  the  APEX  method  plays  a useful  role.  This  is  the  marching  backwards 
in  time  option  characteristic  of  class  G psfs,  which  allows  for  deconvolution  to  be 
performed  in  slow  motion.  Robustness  is  a third  important  property  of  the  APEX 
method,  allowing  detection  of  multiple  psfs  capable  of  significant  sharpening.  This 
substantially  increases  the  probabilities  of  finding  a useful  candidate. 

On  the  theoretical  side,  this  paper  raises  new  questions.  The  first  of  these  is 
the  existence  of  several  useful  psfs.  as  demonstrated  in  the  Sydney  image  in  Figure 
I.  This  phenomenon  warrants  further  investigation.  A second  question  concerns 
the  important  role  Levy  psfs  appear  to  play  in  numerous  imaging  systems.  The 
discussion  in  section  2 has  surveyed  inferences  of  stable  laws  that  have  been  made 
from  mtf  measurements.  Development  of  methods  of  analyzing  imaging  systems  that 
can  rigorously  establish  such  laws,  and  predict  the  Levy  exponent  ft,  would  be  a major 
advance. 

Reconciling  the  results  of  section  2 with  the  behavior  of  large  classes  of  images 
raises  additional  questions.  Electronic  imaging  psfs  h(x,y)  are  found  to  have  expo- 
nents ft  > 0.5  in  most  cases,  so  that  log /;,(£,  0)  = — oj^|2'^  is  a monotone  decreasing 
concave  function  on  f > 0.  However,  as  illustrated  in  Figure  3,  all  images  g{x.  g)  used 
in  this  paper  are  such  that  global  behavior  in  log  | <?*(£,  0)|  is  generally  monotone 
decreasing  and  convex.  Another  large  class  of  images  with  this  convexity  property, 
the  class  W,  was  exhibited  in  [9].  When  such  images  are  APEX-fitted  with  a Levy 
psf  in  the  manner  shown  in  Figure  5(B),  a value  of  ft  < 0.5  is  inevitably  detected. 
An  average  value  of  ft  = 0.23  was  found  for  the  six  images  in  Figures  4,  C,  11.  15, 
16,  and  17,  and  significantly  lower  values  were  found  for  the  remaining  three  images 
in  Figures  12,  13,  and  14.  A possible  partial  explanation  for  this  discrepancy  is  pro- 
vided by  the  Sydney  experiment  in  Figure  1.  There,  the  APEX  method  detected 
several  useful  psfs  with  values  of  ft  smaller  than  the  value  that  was  used  to  blur  the 
image.  The  detected  /i- values  in  the  above  nine  images  may  likewise  underestimate 
the  true  imaging  system  /Lvalues.  An  entirely  different  scenario  may  be  that  the 
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APEX  met  hod  provides  generic  low  exponent  Levy  psfs  capable  of  enhancing  a wide 
variety  of  images,  independently  of  the  imaging  physics  that  created  them.  Other 
generic  enhancement  techniques  have  been  used  for  some  time  in  image  processing, 
[28,  Chap.  10].  More  recent  approaches  based  on  nonlinear  diffusion  equations  are 
also  intended  as  generic  enhancement  methods,  [10].  However,  such  methods  require 
large  numbers  of  iterations  and  are  not  well  suited  for  real-time  processing  of  complex 
high  resolution  imagery. 

Whatever  may  be  the  reasons  behind  it,  the  effectiveness  of  the  APEX  method 
on  many  types  of  images  is  undeniable,  and  the  method  is  a useful  addition  to  the 
image  processing  toolbox. 
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